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^ ! Abstract 

\ We investigate the general relativistic collapse of spherically symmetric, massless spin-^ fields 

at the threshold of black hole formation. A spherically symmetric system is constructed from 
^H^. two spin-^ fields by forming a spin singlet with no net spin-angular momentum. We study the 

system numerically and find strong evidence for a Type II critical solution at the threshold between 
^ ' dispersal and black hole formation, with an associated mass scaling exponent 7 ~ 0.26. Although 



the critical solution is characterized by a continuously self-similar (CSS) geometry, the matter 
fields exhibit discrete self-similarity with an echoing exponent A ~ 1.34. We then adopt a CSS 
ansatz and reduce the equations of motion to a set of ODEs. We find a solution of the ODEs that 
is analytic throughout the solution domain, and show that it corresponds to the critical solution 
found via dynamical evolutions. 
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I. INTRODUCTION 



Beginning with an investigation of the spherically symmetric collapse of a massless scalar 
field Pj, many studies of gravitational collapse have established that the threshold of black 
hole formation is characterized by critical phenomena analogous to those that accompany 
phase transitions in statistical mechanical systems. Briefly, (for a detailed review, see 
black-hole threshold phenomena arise from the consideration of families of spacetimes, (gen- 
erally containing one or more matter fields), which are labeled by a family parameter, p, that 
controls the degree of self-gravitation in the spacetime. These families typically describe the 
implosion of an initially in-going concentration of matter-energy. For small values of p, the 
energy implodes in an essentially linear fashion, re-emerges and disperses to large distances. 
In contrast, for large p, the implosion results in black hole formation, with some fraction 
of the initial mass of the system trapped within a horizon. The threshold of black hole 
formation is defined by the specific (critical) parameter value, p*, at which a black hole first 
makes an appearance. 

Empirically, and quite generically, a number of intriguing features are seen in the near- 
critical regime. These include the emergence of specific critical solutions with scaling, or 
time-translational, symmetries, scaling laws of dimensionful quantities such as the black-hole 
mass, and universality in the sense that these features do not depend on the details of the 
family used to generate the critical solution. There is now a relatively complete, though 
non-rigorous, understanding of this phenomenology. First, for a given matter model and 
symmetry restrictions (spherical, or axial symmetry for example), the black hole threshold 
apparently defines specific solutions of the coupled matter-Einstein equations. These solu- 
tions are essentially unique, up to certain rescalings, or at least are isolated in the overall 
solution space of the model. Second, these critical solutions, although unstable essentially by 
construction, tend to be minimally unstable in the sense of possessing a single (or perhaps a 
few) unstable modes in perturbation theory. The Lyapunov exponents associated with these 
modes can be immediately related to the exponents determined from empirically measured 
scaling laws. 

In this paper we study the critical collapse of a massless spin-i Dirac field coupled to 
the general relativistic gravitational field, within the context of spherical symmetry. Since a 
single spinor field cannot be spherically symmetric, we consider the incoherent sum of two 
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independent fields so that the superposition has no net spin-angular momentum. Through 
direct solution of the 1+1 dimensional PDEs governing the dynamics, we demonstrate the 
existence of so-called Type II behavior in the model, in which black hole formation turns on 
at infinitesimal mass, and the critical solution is self-similar. In the current instance, the self- 
similarity is somewhat novel in that individual components of the Dirac field are discretely 
self-similar, but the overall geometry is continuously self-similar. As expected for this type 
of critical solution, we find a black hole mass-scaling law for solutions in the super-critical 
regime. We then directly construct the threshold solution using an appropriate self-similar 
ansatz for the geometry and matter fields and demonstrate, good agreement between it and 
the PDE solution. 



II. FORMALISM 

We consider a spherically symmetric system of spin-| fields. We use the ADM formalism 
see for details), adopt units in which G = c = h = 1, and express the metric in polar-areal 
coordinates: 

ds'^ = -a{t, rfdt^ + a(t, rfdr"^ + r'^dO'^ + sin^ Odcj)'^ . (1) 

The coordinate r is the areal radius defined as (^4/47?)^, where A is the proper area of a 
constant-r 2-sphere. The functions r) and a{t,r) are to be determined using a subset of 
the 3 + 1 form of Einstein's equations, as described in more detail in Section III (^1 

Before discussing how we separate the radial and angular dependences in our system to 
brm a spin singlet, we begin with a brief review of spinors in curved spacetime (see 0|, 
^, and p]). The evolution of a massless, spin-| field coupled to gravity is governed by the 
curved space Dirac equation: 

j'^V^ij = 0. (2) 

The curved space 7-matrices, 'j^ satisfy 

9'"'t = l{r,Y} (3) 

where U is the 4x4 identity matrix and g^'^ is the inverse metric. In flat spacetime, we 
have: 
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A particular choice of 7° that satisfy Q with our metric signature (— , +, +, +) is: 

1 \ . / (7 A 

7° = M ' 7^ = M ■ (5) 

\^ -1 y y-a^ J 

Here, the index j ranges over the spatial values 1, 2, 3 and the are the Pauli spin matrices, 
namely 

-(;:)•-(::;)■-(::.) 

The general 7-matrices are related to their flat, Cartesian counterparts by 

r = v/T (7) 

where there is an implied summation over the values 0, 1,2,3 of the "fiat", Latin index a, 
and the VJ^ are known as vierbein. 

The derivative operator in equation (j2I) is a spinor covariant derivative with spinor affine 
connections, F^. It acts in the following way on spinors, 

and on 7-matrices, 

^,1" = ^7^^ + r%A7' - + 7^r,. (9) 

However, it reduces to the usual covariant derivative when acting on tensors. We choose the 
spinor connections so that 

V^7'^ = 0. (10) 
It can then be shown that the T ^ take the form 

rM = -^[r,7']K'^v^H.. (11) 

We also note that when taking the covariant derivative of the vierbein, VJ^, only one Christof- 
fel connection appears, since there is only one curved, tensor (Greek) index. 



A. Representation 



Having fixed the form of the spherically symmetric metric (^, we are ready to find a set 
of 7-matrices that satisfy equation (jHl). We choose as our representation: 



The spinor connections are 



1 sin 6 



7°f, 



la' 

1 a 
2a 



0-77', 
2 a 



— 7^7"^ H — cos6'7^7"^. 
a 2 



(13) 



where dots and primes denote differentiation with respect to t and r respectively. We note 
that we have complete freedom to choose any set of 7^^ we wish provided they satisfy equation 
(jni), and that our specific choice is made so that the Dirac equation can be easily separated 
into radial and angular parts. 

Before proceeding to this separation, we introduce a further simplification based on the 
fact that we are dealing with a massless spin-i field. Mathematically, such a field has a 
particular chirality (circular polarization); we adopt left-handed chirality which is expressed 
as: 

(1 - Z7^) ^p = 0. (14) 

Here, 7^ is defined by 

7^ = ffff. (15) 
Equation (fT^ can be satisfied by taking 

Substitution of this form for the spinor into equation Q yields two identical sets of equations, 
each coupling the spinor components, ipi and ip2- We of course only need to solve one set of 
equations for these variables, so we are left with two equations instead of the original four. 
We now perform a separation of variables on the spinor components by writing: 



(16) 



F{t,r)Hi( 



(17) 



With this new choice of variables, and with our previously chosen representation of the 
7^ (fT^ . the Dirac equation separates into a part that depends on [t, r) (which we will refer 
to as the "radial" part), and a part that depends on (6*, 0): 

ir ir a' / -F/G \ ir / F'/G \ ir a' / F/G \ 

" \^ G/F j y Qjp j ay q> j p j 2 aa ^ _Q/p j 

We note that although the factor {r^/a)~^ in ()17|) is not necessary for the separation of 
variables, it simplifies matters by reducing the number of terms in the resulting equation. 
Of particular importance is the elimination of a time derivative of a(t, r) that would make 
numerical solution of the resulting system of PDEs somewhat more involved. 

Considering our separated equation, we observe that since any change m. 9 oi (j) cannot 
change the value of the (t, r) part of ()18j) . the {9, 0) part must be a constant. At this point, if 
our goal was to simply remove the angular dependence from the Dirac equation, we would be 
done. By replacing the angular part of (fTHj) by a constant we would be restricting ourselves 
to some spinor that is an eigenfunction of the angular operators in the Dirac equation. In 
fact only one of its eigenvalues, rather than the precise form of the angular eigenfunction, 
would need to be known. However, our goal is not only to eliminate the angular dependence 
of our equation of motion, but also to have a matter source that generates a spherically 
symmetric spacetime. An individual spinor is not a spherically symmetric object (it always 
has a spin-angular momentum that breaks this symmetry) and therefore cannot by itself 
produce such a spacetime. What we require are multiple spinors where all the individual 
spin-angular momenta counterbalance each other so the system has no net spin. We will 
use two spinors, for simplicity, but any even number of the appropriate spinors could also 
be used. The sphedeal.y sym,„etnc stress-energy tensor for the system. r„„. .s found fro,n 
the sum of the stress-energy tensors of the individual spinor fields |7|]. 

Evaluating the right hand side of equation p9|) does require the angular eigenfunctions that 
we will now compute. 
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B. Equations of Motion 



Setting the angular part of equation (fTHj) equal to a constant, n, gives: 

1 



i 


d 


d 


sin6' 


d(j) 


~ m 


i 


d 


d 


sin 6 


d(j) 





■ cote 



1 

2 



H2 = -nHi, 
Hi = nH2, 



(20) 
(21) 



where we have multiplied the first and second components of the angular terms in (fTTj) by 
—Hi and H2 respectively, so that the bracketed terms in the above expression are the raising 
and lowering operators, S (eth) and 8 (ethbar), respectively. These operators act on the 
spin weighted spherical harmonics, ^im, (see js| and j^) in the following way: 



diXlm) = v/(/-s)(/ + S+l)(,+il^J 



BiYim) = -V(^ + S)(/-S + I)(,_il^j. 



Our functions Hi and H2 have the spin weights s = ±i: 



Hii 
H2U 



(22) 
(23) 

(24) 
(25) 



To form a spin singlet, all we require is one spinor constructed from lYi 1, _iyii, and 

222 222 

one spinor from \Yi_i, _iyi„i. These spin weighted spherical harmonics are: 
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2 22 



1^11 

2 22 
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-^e^<^/^in-, 
2' 

-^e^^^/^cos-, 
2' 

1 ,e-^<^/2^osl 



^ e-^'t'l^siJ- 



'2ti 2 

and are solutions to ()20p and (j?H) for n = — 1. Thus, for our two spinor fields, we have: 



(26) 



,20/2 



( F{t,r) sm{9/2) \ 
G{t,r) cos{e/2) 
F{t,r) sm{9/2) 
\G{t,r) cos(^/2) / 



(27) 



-i(t>/2 



( F{t,r)cos{e/2) \ 
-G{t,r) sm{e/2) 
F{t,r) cos(^/2) 
\-G{t,r) sm{9/2)J 



(28) 



From the above results, we can now derive the following radial equations of motion from ()18p: 



G2 



a 



a 



Go 



air 



-\ -dr 



a 



a \\ a I r 



(29) 

(30) 
(31) 
(32) 



a \\ a J r 

where we have written the complex functions, F{t, r) and G{t, r) in terms of real functions 
via F(t,r) = Fi{t,r) + iF2{t,r) and G{t,r) = Gi(t,r) +iG2it,r). 



C. Geometry 

Although the spinors ()27p and (PHjl both yield the same radial equations of motion ()29|) - 
dnil), they have different stress-energy tensors. We calculate the stress-tensor for each field 
individually using: 

T^,u = [tpl(f,^u)'4' - ^{iM^) lv)^\ (33) 
where the Dirac adjoint of t/', t/' is defined by 

Here, A is the so-called Hermitizing matrix, needed in the computation of real-valued ex- 
pressions (such as the current density or, in this case, the stress-energy tensor) from the 
complex-valued spinors. It is to be chosen so that both A and iA^^ are Hermitian, and we 
take A = —ij^. 

Computing T^^, for each of the two spinors (j77j) and (PS|) and summing the results yields 
the following non- vanishing components of the spherically symmetric stress-energy tensor: 

Ttt = ("^^"^^ ~ + G1G2 — GiG-i^ 
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a 



F1F2 — F1F2 + G1G2 — G1G2 H — {F1F2 ~ F1F2 + G'iG2 — G1G2 



a 

1 

T00 = {FiGi + F2G2) 



2) 



-^^^ ~ o 2 (-^1-^2 ~ -^1-^2 + G'iG2 — GiG'2) 



sin^ B 

Th=o {FiGi + F2G2). (34) 

Contracting the stress-energy tensor gives 

T/ = 0, (35) 

which is expected since the massless Dirac system is conformally invariant. Having computed 
a stress-energy tensor that will generate a spherically symmetric spacetime, we can now write 
down the Einstein equations that will fix a{t, r) and a{t, r). 

Due to our choice of coordinates, the Hamiltonian constraint and the slicing condition 
comprise a sufficient subset of the Einstein equations to be used in our numerical solution. 
The Hamiltonian constraint is 

- + = ^ (2aFiGi + 2aF2G2 + rF^F^ - rF[F2 + rG\G2 - rGiG'2) (36) 

and is treated as an equation for a{t,r). We note that the momentum constraint, 

20- 

d = — {F;F2 - F1F2' + G[G2 - GiG'2) . (37) 

also yields an equation for a (an evolution equation) that we use as a means to check the 
consistency of our equations, both at the analytic level, and during numerical evolutions. 
We note that in both equations ()36|) and ()37|) . time derivatives of F's and G's have been 
eliminated using the equations of motion (j29|l - (|32|) . 

The slicing condition, which fixes a{t,r), is derived from the evolution equation for K^g 
and the fact that for polar slicing we have 

K = K\ = K\ + 2K\ = K\ 

since K^g{t, r) = 0. To maintain K^eit, r) = for all time, we impose K^g{t, r) = 0, which 
yields 

o' - 1 2 

^ = - (F1F2' - F[F2 + G[G2 - G.G'^) . (38) 

a 2r r 
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III. NUMERICS AND RESULTS 



F^{t,0) 


= 


F2{t,0) 


= 




= 


G2{t,0) 


= 



The above equations of motion were solved using a Crank-Nicholson update scheme, 
standard 0{h'^) spatial derivatives, and Berger-Oliger style adaptive mesh refinement (see 
[l^ ) on the computational domain, < r < rmax, t > 0. To achieve stability, high frequency 
modes were damped using Kreiss-Oliger dissipation [ll|. At r = 0, the following regularity 
conditions were enforced: 



(39) 



At r = rjnax, outgoing wave conditions are imposed on the matter fields: 

im = -drFu 

dtGi — —drGi, 
9tG2 = —drG2- 

We also have 

a(t,0) = l, (41) 

which follows from the demand of regularity (local fiatness) at r = 0. At each time step, the 
Hamiltonian constraint is integrated outwards from r = using a point-wise Newton method. 
The shcing condition (|nHj) is solved subject to the outer boundary condition a(t,rmax) = 
l/a(t,rmax), so that, provided that all of the matter remains within the computational 
domain, coordinate time and proper time coincide at spatial infinity. We note that the 
more natural normalization choice, from the point of view of critical collapse, is a{t, 0) = 1. 
However, computationally, this choice would vitiate our ability to compute with a fixed 
Courant factor, At/Ar, particularly in the near-critical regime, and would thus unnecessarily 
complicate the numerical solution. 

The results presented below were generated from the study of three distinct parametrized 
families of initial data, which we refer to as Gaussians, spatial derivatives of Gaussians, and 
kink-anti-kink. In all cases, the initial datasets were such that the pulses were essentially 
ingoing with positive energy at t = 0. The specific forms used are as follows: 
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Gaussian: 

Fi = 
F2 = 

Derivative of Gaussian: 

Fi = 
F2 = 

G^=p ((ro - r)/(252))e-('--«)V4<5^ 
G2=p ((ro - r - (5)/(252))e-('-''o+^)'/4^' 

Kink anti-Kink: 

Fi = 
F2 = 

= (p/2)(tanh((r - ro)/6) - tanh((r - 2ro)/5)) 
G2 = (p/2)(tanh((r - tq + - tanh((r - 2ro + 6)/ 6)) 

In each case the family parameter, p, can be used to control whether the mass-energy of 
the system collapses to form a black hole, or if it implodes through the center (r = 0) and 
disperses to spatial infinity. As discussed in the introduction, and in a process now familiar 
from many other studies of critical collapse, as p is tuned to the black hole threshold p = p*, 
the single unstable mode associated with the critical solution is "tuned away" to reveal the 
critical solution per se. 

In the current case, we find strong indications from such studies that the critical solution 
for the spherically symmetric EMD (Einstein-massless-Dirac) system describes a continu- 
ously self-similar (CSS) geometry. Typical evidence for this claim is shown in Fig. El which 
displays near-critical evolution of the scale invariant quantity, a. In contrast, the Dirac fields 
themselves appear to be discretely self-similar (DSS) (see Fig. |^. 

As usual, associated with the Type II critical solution is a scaling law for the black hole 
mass near criticality: 

M,h^\p-p*\\ (42) 
11 



and, as expected, we find strong evidence that the scaling exponent, 7, is universal in that it 
is independent of the family of initial data used to generate the critical solution. The values 
of 7 computed for the different families are summarized in Table HI and we note that there 
is uncertainty in the third digit of the quoted values. 

TABLE I: Scaling exponent, 7, associated with the three famihes of initial data described in the 
text. 



Family 


7 


Gaussian 


0.258 


derivative of Gaussian 


0.259 


kink anti-kink 


0.257 



The data in Fig. C] does not oscillate about the fit line since the geometry is continuously 
self-similar. This lack of regular oscillations about the fit line is shown more clearly in Fig. |21 



IV. SELF-SIMILAR ANSATZ 



Given the numerical evidence suggesting the existence of a self-similar solution at the 
black hole threshold, we proceed to an ab initio computation of the critical solution based 
on the application of a self-similar ansatz to our system. The development here closely 



12 1, who considered the case of a 



parallels the work done by Hirschmann and Eardley 
massless, complex scalar field. 

By definition, a self-similar spacetime has a homothetic Killing vector, ^, that obeys 

^agf^,^ = '^g,j^u (43) 

where the factor 2 is simply a matter of convention. We want to define coordinates, (r, x) 
that are adapted to this self-similar symmetry. Specifically with r adapted to the vector 
field (pSj) can be written as: 

d 

x) = 2g^^{r, x). (44) 

Performing a separation of variables on the metric tensor and then solving (j44p for the 
r-dependent part yields 

9iiu{r,x) = e'^''g^^{x), (45) 
12 



-30 -20 -10 

ln|p — p* 



r = In 



FIG. 1: Plot of the observed mass scaling near critical! ty for the case of the Gaussian family. The 
measured scaling exponent is 7 = 0.258, with uncertainty in the third digit. As we tune arbitrarily 
close to the critical parameter, p*, black holes of arbitrarily small masses are formed, indicative of 
a Type II critical solution. The small irregularities visible in the plot are shown and discussed in 
more detail in Fig. |2l 

where gp.v{x) is the part of the metric that depends only on x. The original coordinates 
(t, r) are related to (r, x) by: 

'f'-t r 

—— , X = . (46) 

L t*-t ^ ' 

The equations of interest will take the same form for any value of the constant factor, 
L, and without loss of generality, we subsequently take L = 1. The time is the value 
of coordinate time to which the self-similar solution asymptotes as it propagates down to 
arbitrarily small scales, and, from the point of view of critical evolution, is a natural temporal 
origin for use in describing the dynamics. Further, as mentioned previously, in analyzing self- 
similar critical solutions, it is natural to adopt a parameterization of the t=constant surfaces, 
such that if: coincides with central proper time. We thus adopt such a normalization, and 
additionally adjust t so that = 0. 
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FIG. 2: Plot of the residuals of the data shown in Fig. ^ with respect to the computed linear 
fit. The absence of regular oscillations indicates that the geometry is not discretely self-similar. 
There are, however, features in the plot, notably the "spikes", that can be explained as follows. 
The matter fields have a discretely self-similar nature, but combine to produce a continuously self- 
similar geometry. Truncation error effects, combined with the fact that our determination of black 
hole mass is not precise, result in a residual imprint of the DSS nature of the Dirac fields in the 
plot. In fact the spikes in this plot occur as the amplitude of the Dirac fields reach an extremum 
near the black hole, periodically in In |p — 

In these new coordinates, the metric (0) becomes 

ds^ = e^^ [{-a{xf + x^a{xf)dT'^ + 2xa{xfdTdx + a{xfdx^ + x^dO"^ + x^ sin^ Odcf)'^] . (47) 

We note that the r coordinate is timehke, and, as can be verified by comparing the right hand 
sides of (HHjl and (jUj), that the functions, a and a are functions of the spacelike coordinate, 
X, alone. 
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In (r, x) coordinates, the spinors ()27|) and ()28|) become 

/ F(r,x) sin(^/2) \ 
g-r gi0/2 G(r,x)cos(^/2) 
2v^xv/aOr) F(r, x) sin(e/2) 

\G{t,x) cos(0/2) / 

/ F(r, a;) cos(^/2) \ 
e-" -^(r, s) sin(^/2) 

2Axv/a(x) F(r,x)cos(^/2) 

\-G{t,x) sm{e/2) J 



(48) 



(49) 



as 



where these expressions were found by transforming the (t, r) parts of (j^7|) and 
scalars. 

In order to find spinor components that are only functions of x, we require knowledge of 
the T dependence of our field quantities. This is determined by performing the coordinate 
transformations on the equations of motion (P^ - ljH^ and the geometric equations ()36|) and 
(IHHj) . and then ascertaining what r dependence is needed in F and G to produce a set of r 
independent ODEs. A suitable ansatz is 

F{t,x) = e^/^e'^^^xiPiix) +iP2{x)) (50) 

G{t,x) = e^/^e'^^x{Q,{x)+tQ2{x)). (51) 

where the exp{iujt) terms reflect the fact that, as the PDE solutions have revealed, we expect 
the matter fields to exhibit discrete self-similarity. Note that u as defined here corresponds 
to 27r/A (see ^^) where A is the echoing exponent originally defined in p,]. Additionally, 
the extra factor of x is introduced to cast the resulting equations in a more convenient form. 
Inserting this ansatz into (jHEI) and (jHH|) . we find 



Pi 



X + a/a 
1 



la / + 1 
2 a 



X 



a 



P' = 

^ X + a/a 



Q'l 
Q2 



■ 1^ ^ la fa^ + l 
2 2 a \ X 



X — a/a [_ 2 
1 



2 a \ X 



X — a/a 



1 la f a^ + 

2 2 a \ X 



Pi + 2aPi(PiQi + P2Q2) + -Q2 

X 



a 

P2 + 2aP2(PiQi + P2Q2) - -Qi 

X 

Qi-2aQ,{P,Q, + P2Q2) + -P2 

X 



a 

Q2 - 2aQ2{PiQi + P2Q2) - -Pi 

X 



(52) 

(53) 
(54) 
(55) 
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^ = + 2x (PiP^ - P[P2 + Q;Q2 - Q1Q2) + 4a (PiQi + P2Q2) (56) 

a' - 1 

- = + 2x (PiP^ - Pi'P2 + g;g2 - QiQ'2) ■ (57) 

We note that we can cast this system into a canonical form suitable for numerical integration, 
by using ()52|) - (jKK|) to eliminate the derivatives of P and Q that appear in the right hand 
side of the equations for a fjKH|) and a (jKTjl . 

V. NUMERICAL SOLUTION OF THE SELF-SIMILAR ODES 

Having rewritten our equations in a coordinate system adapted to self-similar symmetry, 
we can now solve the resulting ODEs to determine what we expect will be the CSS critical 
solution seen at the black hole threshold in the EMD model. Following Hirschmann and 



Eardley [1^, we use a multi-parameter shooting method to integrate the equations subject 
to regularity and analyticity conditions. 

We first observe that the system (|H^ - (jH7j) has singularities at x = and at x = X2 = a/a 
(the similarity horizon). Of the infinitely many solutions to the ODEs, we seek one that is 
analytic at both of these points, as the CSS solution found via solution of the PDEs has this 
property. Our unknown problem parameters ( "shooting" parameters) include the values of 
some of the fields at the origin, the value of u appearing in the ansatz (fSn|) - (f^ . as well as 
the value of X2 (the position of the similarity horizon). Following Eardley and Hirschmann 
we shoot outwards from x = and inwards from X2 = a /a, comparing solutions at some 
intermediate point xi. This process is automated by starting from some initial guess, then 
using Newton's method to determine the shooting parameters for subsequent iterations. In 
the Newton iteration, we use the square of the differences of the values of the functions 
and their derivatives at xi (computed from the inwards and outwards integrations) as the 
goodness-of-fit indicator, which is driven to as the iteration converges. 

At X = we have the following: 

Pi(0)=0 
A(0) = -Qo 
Qi(0) = Qo 
Q2(0) = 
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a(0) = 1 
a(0) = 1. 

Regularity at the origin gives Pi = Q2 and P2 = —Qi- We use the global U{1) invariance of 
our system ()52j) - (jH7j) to set Pi = 0. This leaves Qo as a shooting parameter. 

As noted previously, the location of the similarity horizon (outer boundary of the inte- 
gration domain) X2 is itself a shooting parameter, and is the value x where a{x)/a{x) = x. 
In the limit x ^ X2 we have the following: 

Piix2) = ^ {-AaQluj - 4:aQlQ2UJ + 2aujQi + Qsa^) 

^2(2:2) = ^ {4:aQiQlu + iaQlu + 2aujQ2 - Qia^) 
Qi{x2) = Qi 
Q2ix2) = Q2 
a{x2) = X2a 
a{x2) = a. 

The shooting parameters at the outer boundary are: X2, Qi{x2), Q2{x2), and a{x2). The 
final shooting parameter is the frequency, u, for a total of six undetermined parameters. 
We find an approximate solution given by 

X2 = 5.6740230 ± 0.0000004 

= 4.698839 ±0.000001 

Qi(0) = 0.747912623 ± 0.000000006 

Qi{x2) = 0.00151341532 ± 0.00000000007 

Q2{x2) = 0.01103266083 ± 0.00000000005 

a(x2) = 1.1183631604 ±0.0000000009. 

where the quoted uncertainty was estimated by solving the system for many different values 
of Xi and observing the changes in the shooting parameters. 
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VI. COMPARISON OF PDE/ODE SOLUTIONS 

In this section we compare the solution computed from the self-similar ansatz, as just 
described, to the near-critical solutions calculated from the full PDEs in the {t, r) coordinate 
system. The ODE solution is the theoretically predicted self-similar solution while the PDE 
solution can be thought of as collected data. For this comparison, we used data from the 
Gaussian family. The idea is to treat the ODE solution as the model function and fit it to the 
PDE data. We perform the fit "simultaneously" at all times by working with the functions 
as 2-dimensional solution surfaces in t and r. This process is automated by starting from 
some initial guess for the fitting parameters, then using Newton's method to determine these 
parameters for subsequent iterations. In the Newton iteration, we use the least squares of 
the two solutions: 

N 

Y^{ur--uf--Y (58) 

as the goodness-of-fit indicator, which is driven towards as the iteration converges. When 
this happens, we compute the Z2-norm of the difference of the two solutions: 

iHb=(^E("r'=-«r'=)j (59) 

as an error estimate. In computing both the least squares and the /2-norm, the solution to 
the ODEs and the solution to the PDEs are treated as A^-element vectors where N is the 
total number of grid points on the 2-dimensional grid in t and r. 

For the case of a{t, r), we use t* as the fitting parameter. Once t* is found, the /2-norm of 
the difference of the two solutions is 0.00159. We display the results of this fitting process 
as a sequence of snapshots in time in Fig. IHl 

Comparing the Dirac fields is a little more involved since there are a number of unspecified 
parameters and phases that must be determined. From the ansatz (j3m) -()51 |) we have: 

Fi(t,r) = , [Piit,r) cos(culn(t* -t) + (/)i) -P2(t,r)sin(cjln(t* -t) + 0i)] 

{t* — t) 2 

F2(t,r) = , [Pi{t,r) sm{uj\n{f - t) + (j)2) + P2{t,r) cos{uj \n{f - t) + (f)2)] 

{t* — t) 2 

Gi(t, r) = 1 [Qi{t, r) cos(o; ln(r - t) + ^2) - <32(t, r) sm{cu ln(r - t) + ^2)] 

(t* — t) 2 



18 



-1-3" 



1.5 



t = 



t = 3.7242 




t = 3.7520 



^ * ODE A 




t = 3.3683 




t = 3,7446 




L-3.7525 




t = 3,6503 




t = 3.7504 




L-3,7527 




In r 



FIG. 3: Evolution of the metric variable, a{t,r), for a slightly supercritical evolution, overlayed 
with the solution to the ODEs. The frames are output logarithmically in central proper time, as 
measured from t*. The function's peak reaches a value of approximately 1.6 and remains there as 
the solution continuously repeats itself on ever smaller scales. 



G2{t,r) 



rAi 



- [gi(t,r)sin(u;ln(t^ -t) + 0i) + g2(t,r)cos(cjln(t* -t) + 0i)] . (60) 



We note that Fi and G2 have the same phase, (pi, while the pair F2 and Gi have the same 
phase 02- This is expected from the coupling of (j^^ -()32j ) . The equations of motion may 
be invariant under changes of these phases, but (jHBj) and (jHHj) are not. In order to have the 
entire system be invariant under changes in the phases, we must have: 

A1A2 ^ 



COS(^</)l - (p2) 

We see that the amplitudes of the fields must change only if the relative phase, 0i — (f>2, 
changes. We merely note this fact for completeness but do not use it to reduce the number 
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FIG. 4: Comparison of Fi{t,r) for a slightly supercritical PDE evolution with the solution found 
from the self-similar ansatz. As in Fig. 01 the frames are output logarithmically in central proper 
time, measured from t*. The function oscillates and discretely repeats itself on smaller and smaller 
scales. The particular combination of Fi(t,r)/ y/r is the natural scaling variable as can be deduced 
by noting that the Dirac fields have units of (Length) and as shown explicitly by the r depen- 
dence of (|5Ujl and (|51|) . The echoing exponent, A, calculated via A = In /to, where u) is computed 
from the ODE solution, is also shown. 

of fit parameters. 

The comparison of the fields as found from the PDEs and ODEs is carried out in much 
the same way as it is done for the metric variable, a{t, r). The goodness-of-fit is again defined 
to be the least squares of the two solutions but this time, the parameter t* is kept fixed and 
the phase (pj and amplitude Aj are used as fitting parameters (j = 1,2). The /2-norm of the 
difference of the solutions for Fi is 0.000195. The /2-norm of the difference of the solutions 
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for Gi is 0.00024. Fig. H] illustrates the results of this comparison of the ODE and PDE 
solutions for Fi. 

VII. CONCLUSIONS 

We have investigated the spherically symmetric Einstein-massless-Dirac system at the 
threshold of black hole formation. We have found strong evidence for a Type II critical so- 
lution, and an associated mass scaling law with a universal exponent 7 ~ 0.26. The solution 
exhibits continuous self-similarity in the geometric variables and discrete self-similarity in 
the components of the Dirac fields, the latter characterized by an echoing exponent A ~ 1.34. 
Using a self-similar ansatz, we then reduced the equations of motion governing the model 
to a set of ODEs whose solution, given appropriate regularity conditions, is in very good 
agreement with the critical solution obtained from the original PDEs. 
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